---
title: "Why do governments tax or subsidize fossil fuels?"
author: "Paasha Mahdavi, Cesar B. Martinez-Alvarez, Michael L. Ross"
date: "October 5, 2021"
output:
  pdf_document: default
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, results = "asis", warning = FALSE, message = FALSE, error = FALSE)
```

```{r}
rm(list = ls())

# Install required packages
pkg <- c("estimatr","stargazer","tidyverse","AER","readstata13","countrycode")
inst <- pkg %in% installed.packages()
if (length(pkg[!inst]) > 0) install.packages(pkg[!inst])
rm(pkg, inst)

library(estimatr)
library(stargazer)
library(tidyverse)
library(AER)
library(readstata13)
library(countrycode)

# Set Working Directory to wherever files are downloaded
# setwd("")
```


```{r}
#Read cross-national data

cross_national_ffs_final <- read.csv("cross_national_ffs_final.csv") %>%
  dplyr::select(-1)

#Read annual data

annualpanel_ffs_final <- read.csv("annualpanel_ffs_final.csv") %>%
  dplyr::select(-1)

#Read monthly data

monthlypanel_ffs <- read.csv("monthlypanel_ffs.csv")
```

```{r tab1}

#
#   Models in main text
#



##Table 1: Cross-Section, Basic Specification (p.12)

model1 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat, data = cross_national_ffs_final)

model2 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model3 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + log(average_oilgas_exports_pc+0.001), data = cross_national_ffs_final)

model4 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + meanfuelexports, data = cross_national_ffs_final)

t1names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
             "Value-Added Tax Rate","Fossil Fuel Dependence","log(Oil and Gas Exports PC)",
             "Fossil Fuel Export Dependence")

stargazer(model1, model2, model3, model4,  
          se = starprep(model1, model2, model3, model4,  se_type = "stata"),
          covariate.labels = t1names,
          keep.stat = c("n","rsq","adj.rsq"),
          title = "Cross-Section / Basic Specification",
          header=FALSE, type='latex')
```

\newpage

```{r tab2}
##Table 2: Cross-Section, Additional Controls (p. 13)

model4 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + avg_gov_debt, data = cross_national_ffs_final)

model5 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model6 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model7 <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)

t2names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model4, model5, model6, model7, 
          se = starprep(model4, model5, model6, model7, se_type = "stata"),
          covariate.labels = t2names,
          keep.stat = c("n","rsq","adj.rsq"),
          title = "Cross-Section / Additional Controls",
          header=FALSE, type='latex')
```

\newpage

```{r tab3}
##Table 3: Instrumental Variables Approach (p. 15)

model1_instrument <- ivreg(meanbmgap2015adj ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2)  + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model2_instrument <- ivreg(meanbmgap2015adj ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model3_instrument <- ivreg(meanbmgap2015adj ~ fuel_income_dependence + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model4_instrument <- ivreg(meanbmgap2015adj ~ fuel_income_dependence + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + + europeanunion + Latitude + landlocked + region2 + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

rob_se <- list(sqrt(diag(vcovHC(model1_instrument, type = "HC1"))),
               sqrt(diag(vcovHC(model2_instrument, type = "HC1"))),
               sqrt(diag(vcovHC(model3_instrument, type = "HC1"))),
               sqrt(diag(vcovHC(model4_instrument, type = "HC1"))))

t3names <- c("Fossil Fuel Dependence",
             "log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Value-Added Tax Rate",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model1_instrument, model2_instrument, model3_instrument, model4_instrument,
          se = rob_se,
          covariate.labels = t3names,
          title = "Instrumental Variables Approach",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex')
```

\newpage

```{r tab4}
##Table 4: Cross-Section, Time-Series: Annual Panel (p. 17)
model_annual1 <- lm(meanbmgap2015adj ~ VAT.Rate + log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

t4names <- c("Value-Added Tax Rate",
             "log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","Fossil Fuel Dependence * Autocracy")

stargazer(model_annual1,  
          se = starprep(model_annual1, se_type = "stata"),
          omit = "factor*",
          covariate.labels = t4names,
          title = "Cross-section Time-series: Annual Panel",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex')
```

\newpage

```{r tab5}
##Table 5: Cross-Section Time Series: Monthly Panel

model_panel2month <- lm(bmgap2015adj ~ one_quarter_before_elections + two_quarter_before_elections + three_quarter_before_elections + four_quarter_before_elections + one_quarter_after_elections + two_quarter_after_elections + three_quarter_after_elections + four_quarter_after_elections + one_quarter_after_leader + two_quarter_after_leader + three_quarter_after_leader + four_quarter_after_leader + oil_discovery_month + one_quarter_discovery + two_quarter_discovery + three_quarter_discovery + four_quarter_discovery + as.factor(country2),  data = monthlypanel_ffs)

t5names <- c("1 Qr Before Elections", 
             "2 Qr Before Elections", 
             "3 Qr Before Elections",
             "4 Qr Before Elections",
             "1 Qr After Elections",
             "2 Qr After Elections",
             "3 Qr After Elections",
             "4 Qr After Elections",
             "1 Qr After Leader Turnover",
             "2 Qr After Leader Turnover",
             "3 Qr After Leader Turnover",
             "4 Qr After Leader Turnover",
             "Oil Discovery Month",
             "1 Qr After Discovery",
             "2 Qr After Discovery",
             "3 Qr After Discovery",
             "4 Qr After Discovery")

stargazer(model_panel2month, 
          se = starprep(model_panel2month, se_type = "stata"),
          omit = "factor*",
          covariate.labels = t5names,
          title = "Cross-section Time-series: Monthly Panel",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex')
```

\newpage

\setcounter{table}{3}
\renewcommand{\thetable}{S\arabic{table}}

```{r s4}

#
#   Models in Appendix
#



##Table S4: Regime Type and Gas Taxes (p. 13)

model1dem <- lm(meanbmgap2015adj ~ meanpolityiv, data = cross_national_ffs_final)

model2dem <- lm(meanbmgap2015adj ~ autocracy_polity, data = cross_national_ffs_final)

model3dem <- lm(meanbmgap2015adj ~ democracy_polity, data = cross_national_ffs_final)

model4dem <- lm(meanbmgap2015adj ~ meanpolyarchy, data = cross_national_ffs_final)

model5dem <- lm(meanbmgap2015adj ~ democracy, data = cross_national_ffs_final)

ts4names <- c("Polity IV (continuous)", "Autocracy (Polity IV)", "Democracy (Polity IV)",
              "Electoral Democracy (V-DEM)", "Democracy (Boix et al.)")

stargazer(model1dem, model2dem, model3dem, model4dem, model5dem,
          se = starprep(model1dem, model2dem, model3dem, model4dem, model5dem,  se_type = "stata"),
          title = "Bivariate relationships between Regime Type and Fossil Fuel Taxes",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts4names)
```

\newpage

```{r s5}
##Table S5: Cross-sectional models: Alternative  measures  of  regime type (Boix) (p. 14)

model1boix <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model2boix <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence + democracy*fuel_income_dependence + meangoveffect, data = cross_national_ffs_final)

model3boix <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + democracy*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model4boix <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + democracy*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)

ts5names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
              "European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "BMR Democracy","Government Effectiveness",
             "Fossil Fuel Dependence * BMR Democracy")

stargazer(model1boix, model2boix, model3boix, model4boix,
          se = starprep(model1boix, model2boix, model3boix, model4boix, se_type = "stata"),
           title = "Cross-sectional models: Alternative measures of regime type",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts5names)
```

\newpage

```{r s6}
##Table S6: Cross-sectional models: Alternative  measures  of  regime type (V-DEM) (p. 15)

model1vdem <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model2vdem <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence + meanpolyarchy*fuel_income_dependence + meangoveffect, data = cross_national_ffs_final)

model3vdem <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + meanpolyarchy*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model4vdem <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + meanpolyarchy*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)

ts6names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
              "European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "Electoral Democracy","Government Effectiveness",
             "Fossil Fuel Dependence * Electoral Democracy")

stargazer(model1vdem, model2vdem, model3vdem, model4vdem, 
          se = starprep(model1vdem, model2vdem, model3vdem, model4vdem, se_type = "stata"),
           title = "Cross-sectional models: Alternative measures of regime type.",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts6names)
```

\newpage

```{r s7}
##Table S7: Cross-sectional models: Alternative  measures  of  regime type (Polity Continuous) (p. 16)

model1politycont <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model2politycont <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence + meanpolityiv*fuel_income_dependence + meangoveffect, data = cross_national_ffs_final)

model3politycont <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + meanpolityiv*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model4politycont <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + meanpolityiv*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)

ts7names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
              "European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "Polity IV","Government Effectiveness",
             "Fossil Fuel Dependence * Polity IV")

stargazer(model1politycont, model2politycont, model3politycont, model4politycont,
          se = starprep(model1politycont, model2politycont, model3politycont, model4politycont, se_type = "stata"),
           title = "Cross-sectional models: Alternative  measures  of  regime type.",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts7names)
```

\newpage

```{r s8}
##Table S8: Cross-sectional models: Alternative measure of region (p. 17)

model1region <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model2region <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect, data = cross_national_ffs_final)

model3region <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model4region <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region_wb, data = cross_national_ffs_final)

ts8names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
              "European Union","Latitude",
             "Landlocked", "Europe + Central Asia", "Latin America + Caribbean",
             "Middle East + North Africa", "North America",
             "South Asia", "Sub-Saharan Africa",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Fossil Fuel Dependence * Autocracy (Polity IV)")

stargazer(model1region, model2region, model3region, model4region,
          se = starprep(model1region, model2region, model3region, model4region, se_type = "stata"),
           title = "Cross-sectional models: Alternative measure of region",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts8names)
```


\newpage

```{r s9}
##Table S9: Cross-sectional models: National oil companies (p.18)

model1noc <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + NOC, data = cross_national_ffs_final)

model2noc <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + powerful_NOC, data = cross_national_ffs_final)

model3noc <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + NOC_cheon, data = cross_national_ffs_final)

ts9names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Value-Added Tax Rate",
             "Fossil Fuel Dependence","Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt", "National Oil Company (Mahdavi)",
             "Influential National Oil Company (Mahdavi)",
             "National Oil Company (Cheon et al.)",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model1noc, model2noc, model3noc,
          se = starprep(model1noc, model2noc, model3noc, se_type = "stata"),
           title = "Cross-sectional models: National Oil Companies",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts9names)
```

\newpage

```{r s10}
##Table S10: Cross-sectional models: Motorization  rate (p. 19)

model1car <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence + average_motorization_rate, data = cross_national_ffs_final)

model2car <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt +  meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + average_motorization_rate, data = cross_national_ffs_final)

model3car <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + average_motorization_rate, data = cross_national_ffs_final)

model4car <- lm(meanbmgap2015adj ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + average_motorization_rate + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)

ts10names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
               "Central Government Debt","Value-Added Tax Rate",
               "Fossil Fuel Dependence","Autocracy (Polity IV)",
               "Government Effectiveness",
               "Motorization Rate",
               "European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Fossil Fuel Dependence * Polity IV")

stargazer(model1car, model2car, model3car, model4car,
          se = starprep(model1car, model2car, model3car, model4car, se_type = "stata"),
           title = "Cross-sectional models: Motorization rate",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts10names)
```

\newpage

```{r s11}
##Table S11: Instrumental Variables models: Alternative instrument (p. 20)

#Alternative Instrument
#Importing Cassidy data to add original instrument

fuel_instrument <- read.dta13("master_oilwealth_EJ_expanded.dta", nonint.factors = TRUE) %>%
  mutate(country2 = countrycode(Country, 'country.name', 'iso3c',warn = FALSE, nomatch = NA)) %>%
  dplyr::select(Country, code3, country2, ln_convCC_mech_1960_pc, ln_land_area_1960_pc, ln_coastline_1960_pc,
         ln_mountains_1960_pc, ln_good_soil_1960_pc, ln_tropics_1960_pc,
         ln_r_op_avgpc_1966_2008, newregion) %>%
  filter(! Country %in% c("Soviet Union","Bosnia"))

cross_national_ffs_final <- cross_national_ffs_final %>%
  left_join(fuel_instrument, by = "country2")

model1_instrument_rob <- ivreg(meanbmgap2015adj ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + avg_gov_debt + ln_land_area_1960_pc + ln_coastline_1960_pc + ln_mountains_1960_pc + ln_good_soil_1960_pc + ln_tropics_1960_pc + newregion| .- fuel_income_dependence + ln_convCC_mech_1960_pc,  data = cross_national_ffs_final)

model2_instrument_rob <- ivreg(meanbmgap2015adj ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2) +  meanvat + avg_gov_debt + autocracy_polity*fuel_income_dependence + meangoveffect + ln_land_area_1960_pc + ln_coastline_1960_pc + ln_mountains_1960_pc + ln_good_soil_1960_pc + ln_tropics_1960_pc + newregion | .- fuel_income_dependence + ln_convCC_mech_1960_pc,  data = cross_national_ffs_final)

model3_instrument_rob <- ivreg(meanbmgap2015adj ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + avg_gov_debt + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt  + ln_land_area_1960_pc + ln_coastline_1960_pc + ln_mountains_1960_pc + ln_good_soil_1960_pc + ln_tropics_1960_pc | .- fuel_income_dependence + ln_convCC_mech_1960_pc,  data = cross_national_ffs_final)


rob_se <- list(sqrt(diag(vcovHC(model1_instrument_rob, type = "HC1"))),
               sqrt(diag(vcovHC(model2_instrument_rob, type = "HC1"))),
               sqrt(diag(vcovHC(model3_instrument_rob, type = "HC1"))))

ts11names <- c("Fossil Fuel Dependence",
               "log(GNI Per Capita)","log(GNI Per Capita Sq)",
               "Value-Added Tax Rate", "Central Government Debt",
               "Autocracy (Polity IV)",
               "Government Effectiveness",
               "log(Land Area PC)","log(Coastline PC)",
               "log(Mountainous Area PC)", "log(Good Soil PC)",
               "log(Tropical Area PC)", "Europe/Northern America/Oceania",
               "Asia", "Latin America/Caribbean",
               "Fossil Fuel Dependence * Polity IV")

stargazer(model1_instrument_rob, model2_instrument_rob, model3_instrument_rob,
          se = rob_se,
           title = "Instrumental Variables models: Alternative instrument",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts11names)
```

\newpage

```{r s12}
##Table S12: Instrumental Variables models: First-stage results (p. 21)

first.out <- lm(fuel_income_dependence~lnendowment_1960_pc, data = cross_national_ffs_final)

first.out_b <- lm(fuel_income_dependence~ln_convCC_mech_1960_pc, data = cross_national_ffs_final)

ts12names <- c("Oil Endowment PC 1960","Basin Type Area PC 1960")

stargazer(first.out, first.out_b, se = starprep(first.out, first.out_b,
                                                se_type = "stata"),
           title = "Instrumental Variables models: First stage results",
          keep.stat = c("n","rsq","adj.rsq","f"),
          header=FALSE, type='latex',
          covariate.labels = ts12names)
```


\newpage

```{r s13}
##Table S13: Cross-sectional Time-series: with and without fixed effects (p. 22)


#No Fixed Effects

model_annual1 <- lm(meanbmgap2015adj ~ VAT.Rate + log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp,  data = annualpanel_ffs_final)

#Fixed Effects

model_annual2 <- lm(meanbmgap2015adj ~ VAT.Rate + log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

ts13names <- c("Value-Added Tax Rate",
             "log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","Fossil Fuel Dependence * Autocracy")


stargazer(model_annual1, model_annual2,
          se = starprep(model_annual1, model_annual2, se_type = "stata"),
          omit = "factor*",
           title = "Cross-sectional Time-series: with and without fixed effects",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts13names)
```

\newpage

```{r s14}
##Table S14: Cross-sectional time-series models:Alternative measures of regime type (p. 23)

#No Fixed Effects

model_annual1_dem <- lm(meanbmgap2015adj ~ VAT.Rate + log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*democracy + goveffect + percent_gdp,  data = annualpanel_ffs_final)


#Fixed Effects

model_annual2_dem <- lm(meanbmgap2015adj ~ VAT.Rate + log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*democracy + goveffect + percent_gdp + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

ts14names <- c("Value-Added Tax Rate",
             "log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Democracy","Government Effectiveness",
             "Central Government Debt","Fossil Fuel Dependence * Democracy")

stargazer(model_annual1_dem, model_annual2_dem,
          se = starprep(model_annual1_dem, model_annual2_dem, se_type = "stata"),
          omit = "factor*",
           title = "Cross-sectional time-series models:Alternative measures of regime type",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts14names)
```

\newpage

```{r s15}
##Table S15: Cross-sectional time-series models: Alternative measures of regime type (p. 24)


#No Fixed Effects

model_annual1_vdem <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*v2x_polyarchy + goveffect + percent_gdp + VAT.Rate,  data = annualpanel_ffs_final)

#Fixed Effects

model_annual2_vdem <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*v2x_polyarchy + goveffect + percent_gdp + VAT.Rate + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)


ts15names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Electoral Democracy","Government Effectiveness",
             "Central Government Debt","Value-Added Tax Rate",
             "Fossil Fuel Dependence * Electoral Democracy")

stargazer(model_annual1_vdem, model_annual2_vdem,
          se = starprep(model_annual1_vdem, model_annual2_vdem, se_type = "stata"),
          omit = "factor*",
           title = "Cross-sectional time-series models: Alternative measures of regime type",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts15names)
```


\newpage

```{r s16}
##Table S16: Cross-sectional time-series models: Alternative measures of regime type (continuous Polity IV) (p. 25)

#No Fixed Effects

model_annual1_politycontin <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*polity2 + goveffect + percent_gdp + VAT.Rate,  data = annualpanel_ffs_final)

#Fixed Effects

model_annual2_politycontin <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*polity2 + goveffect + percent_gdp + VAT.Rate + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

ts16names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Polity IV (Continuous)","Government Effectiveness",
             "Central Government Debt","Value-Added Tax Rate",
             "Fossil Fuel Dependence * Polity IV (Continuous)")

stargazer(model_annual1_politycontin, model_annual2_politycontin,
          se = starprep(model_annual1_politycontin, model_annual2_politycontin, se_type = "stata"),
          omit = "factor*",
           title = "Cross-sectional time-series models: Alternative measures of regime type. ",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts16names)
```

\newpage

```{r s17}
##Table S17: Cross-sectional time-series models: Motorization rate (p. 26)

#No Fixed Effects

model_annual1_car <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + motorization_rate + VAT.Rate,  data = annualpanel_ffs_final)

#Fixed Effects

model_annual2_car <- lm(meanbmgap2015adj ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + motorization_rate + VAT.Rate + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

ts17names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt", "Motorization Rate",
             "Value-Added Tax Rate",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model_annual1_car, model_annual2_car,
          se = starprep(model_annual1_car, model_annual2_car, se_type = "stata"),
          omit = "factor*",
           title = "Cross-sectional time-series models: Motorization rate.",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts17names)
```

\newpage

```{r s18}
##Table S18: Cross-Sectional Time Series (monthly) with and without fixed effects (p. 27)


model_panel2month_fe <- lm(bmgap2015adj ~ one_quarter_before_elections + two_quarter_before_elections + three_quarter_before_elections + four_quarter_before_elections + one_quarter_after_elections + two_quarter_after_elections + three_quarter_after_elections + four_quarter_after_elections + one_quarter_after_leader + two_quarter_after_leader + three_quarter_after_leader + four_quarter_after_leader + oil_discovery_month + one_quarter_discovery + two_quarter_discovery + three_quarter_discovery + four_quarter_discovery + as.factor(country2),  data = monthlypanel_ffs)

model_panel2month <- lm(bmgap2015adj ~ one_quarter_before_elections + two_quarter_before_elections + three_quarter_before_elections + four_quarter_before_elections + one_quarter_after_elections + two_quarter_after_elections + three_quarter_after_elections + four_quarter_after_elections + one_quarter_after_leader + two_quarter_after_leader + three_quarter_after_leader + four_quarter_after_leader + oil_discovery_month + one_quarter_discovery + two_quarter_discovery + three_quarter_discovery + four_quarter_discovery,  data = monthlypanel_ffs)


ts18names <- c("1 Qr Before Elections", 
             "2 Qr Before Elections", 
             "3 Qr Before Elections",
             "4 Qr Before Elections",
             "1 Qr After Elections",
             "2 Qr After Elections",
             "3 Qr After Elections",
             "4 Qr After Elections",
             "1 Qr After Leader Turnover",
             "2 Qr After Leader Turnover",
             "3 Qr After Leader Turnover",
             "4 Qr After Leader Turnover",
             "Oil Discovery Month",
             "1 Qr After Discovery",
             "2 Qr After Discovery",
             "3 Qr After Discovery",
             "4 Qr After Discovery")

stargazer(model_panel2month, model_panel2month_fe,
          se = starprep(model_panel2month, model_panel2month_fe, se_type = "stata"),
          omit = "factor*",
           title = "Cross-Sectional Time Series (Monthly) with and without fixed effects.",
          keep.stat = c("n","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts18names)
```

\newpage

```{r s19}
##Table S19: Diesel analysis: cross-sectional base specification (p. 28)

model1_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat, data = cross_national_ffs_final)

model2_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + fuel_income_dependence, data = cross_national_ffs_final)

model3_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + log(average_oilgas_exports_pc+0.001), data = cross_national_ffs_final)

model4_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + avg_gov_debt + meanvat + meanfuelexports, data = cross_national_ffs_final)


ts19names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)","Central Government Debt",
             "Value-Added Tax Rate","Fossil Fuel Dependence","log(Oil and Gas Exports PC)",
             "Fossil Fuel Export Dependence")

stargazer(model1_diesel, model2_diesel, model3_diesel, model4_diesel,
          se = starprep(model1_diesel, model2_diesel, model3_diesel, model4_diesel,  se_type = "stata"),
           title = "Diesel analysis: cross-sectional base specification",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts19names)
```

\newpage

```{r s20}
##Table S20: Diesel analysis: cross-sectional full specification (p. 29)

model5_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + avg_gov_debt, data = cross_national_ffs_final)

model6_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model7_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt, data = cross_national_ffs_final)

model8_diesel <- lm(mean_bmgap_diesel_spotprice_la ~ log(meangdppcatlas) + I(log(meangdppcatlas)^2) + meanvat + fuel_income_dependence + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + europeanunion + Latitude + landlocked + region2, data = cross_national_ffs_final)


ts20names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Value-Added Tax Rate","Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model5_diesel, model6_diesel, model7_diesel, model8_diesel,
          se = starprep(model5_diesel, model6_diesel, model7_diesel, model8_diesel, se_type = "stata"),
           title = "Diesel analysis: cross-sectional full specification.",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts20names)
```

\newpage

```{r s21}
##Table S21: Diesel analysis: instrumental variable specification (p. 30)

model1_instrument_diesel <- ivreg(mean_bmgap_diesel_spotprice_la ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2)  + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model2_instrument_diesel <- ivreg(mean_bmgap_diesel_spotprice_la ~ fuel_income_dependence  + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model3_instrument_diesel <- ivreg(mean_bmgap_diesel_spotprice_la ~ fuel_income_dependence + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

model4_instrument_diesel <- ivreg(mean_bmgap_diesel_spotprice_la ~ fuel_income_dependence + log(meangdppcatlas) + I(log(meangdppcatlas)^2) + autocracy_polity*fuel_income_dependence + meangoveffect + avg_gov_debt + + europeanunion + Latitude + landlocked + region2 + meanvat | .- fuel_income_dependence + lnendowment_1960_pc,  data = cross_national_ffs_final)

rob_se_diesel <- list(sqrt(diag(vcovHC(model1_instrument_diesel, type = "HC1"))),
               sqrt(diag(vcovHC(model2_instrument_diesel, type = "HC1"))),
               sqrt(diag(vcovHC(model3_instrument_diesel, type = "HC1"))),
               sqrt(diag(vcovHC(model4_instrument_diesel, type = "HC1"))))


ts21names <- c("Fossil Fuel Dependence",
             "log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt","European Union","Latitude",
             "Landlocked", "Asia + Pacific", "Europe + North America", 
             "Former USSR", "Latin America + Caribbean", "Middle East",
             "Value-Added Tax Rate",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model1_instrument_diesel, model2_instrument_diesel, model3_instrument_diesel, model4_instrument_diesel,
          se = rob_se_diesel,
           title = "Diesel analysis: instrumental variable specification. ",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts21names)
```

\newpage

```{r s22}
##Table S22: Diesel analysis: cross-sectional time-series specification (p. 31)

#No Fixed Effects

model_annual1_diesel <- lm(bmgap_diesel_spotprice_la ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + VAT.Rate,  data = annualpanel_ffs_final)

#Fixed Effects

model_annual2_diesel <- lm(bmgap_diesel_spotprice_la ~ log(NY.GNP.PCAP.CD) + I(log(NY.GNP.PCAP.CD)^2) + fuel_income_dependence*autocracy_polity + goveffect + percent_gdp + VAT.Rate + as.factor(country2) + as.factor(year),  data = annualpanel_ffs_final)

ts22names <- c("log(GNI Per Capita)","log(GNI Per Capita Sq)",
             "Fossil Fuel Dependence",
             "Autocracy (Polity IV)","Government Effectiveness",
             "Central Government Debt",
             "Value-Added Tax Rate",
             "Fossil Fuel Dependence * Autocracy")

stargazer(model_annual1_diesel, model_annual2_diesel,
          se = starprep(model_annual1_diesel, model_annual2_diesel, se_type = "stata"),
          omit = "factor*",
           title = "Diesel analysis: cross-sectional time-series specification.",
          keep.stat = c("n","rsq","adj.rsq"),
          header=FALSE, type='latex',
          covariate.labels = ts22names)
```
